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Abstract 

In the present paper we present a finite element approach for option 
pricing in the framework of a well-known stochastic volatility model with 
jumps, the Bates model. In this model the asset log-returns are assumed 
to follow a jump-diffusion model where the jump component consists of a 
Levy process of compound Poisson type, while the volatility behavior is 
described by a stochastic differential equation of CIR type, with a mean- 
reverting drift term and a diffusion component correlated with that of the 
log-returns. Like in all the Levy models, the option pricing problem can 
be formulated in terms of an integro-differential equation: for the Bates 
model the unknown F{S, V, t) (the option price) of the pricing equation 
depends on three independent variables and the differential operator part 
turns out to be of parabolic kind, while the nonlocal integral operator is 
calculated with respect to the Levy measure of the jumps. In this paper 
we will present a variational formulation of the problem suitable for a finite 
element approach. The numerical results obtained for european options 
will be compared with those obtained with different methods. 

* MOX- ModeUistica e Calcolo Scientifico 
Dipartimento di Matematica "F. Brioschi" 
Politecnico di Milano 
via Bonardi 9, 20133 Milano, Italy 
edie . miglioOpolimi . it 
" Dipartimento di Matematica "F. Brioschi" 
Politecnico di Milano 
via Bonardi 9, 20133 Milano, Italy 
carlo . sgarraSpolimi . it 

Keywords: Option Pricing, Stochastic Volatility Models, Levy Processes, Partial- 
Integro-Differential-Equations, Finite Element Methods 

AMS Subject Classification: 62P05, 60G35, 45K05, 65L60 



* Submitted for publication to Computing and Visualization in Science 



1 



1 Introduction 



A huge effort lias been made in tlie last few years in order to overcome the 
intrinsic limitations of the Black-Scholes model. Although it has been a great 
success as a first attempt to provide an evaluation for financial derivatives, it 
was soon clear that it's description of financial market behavior was not satis- 
factory. Very well known observed empirical features of the log prices were not 
correctly described by this model: heavy tails, volatility clustering, aggregational 
gaussianity are some peculiarities that cannot be explained on the basis of the 
lognormal assumption on which the Black-Scholes model stands. The volatility 
smile is another relevant phenomenon that cannot be explained on the basis of 
a Black-Scholes description. Several different approaches have been exploited 
in order to give a more satisfactory description of financial markets, but the 
main contributions in this direction can be grouped in two different classes of 
models, the so called stochastic volatility models and the models with jumps. 
An extended literature is available on both kind of approaches and they both 
give a more realistic description of the prices evolution in financial markets, but 
separately considered they perform very well only in some situations. While 
jump models can succesfully reproduce the volatility smiles on short term ma- 
turity ranges, stochastic volatility models give a better description of the same 
phenomenon on long maturity terms. This has naturally led to the introduc- 
tion of more complicated, but more realistic models in which both features of 
stochastic volatility and jumps can be present. The three more popular models 
in which this integration of jumps and stochastic volatility has been performed 
are the BNS model introduced by O. BarndorfF-Nielsen and N. Shepard [3], [1], 
the model introduced by Bates in [5], and the time-changed Levy models intro- 
duced by Carr, D. Madan, H. Geman and M. Yor in [18| . While in the former 
the volatility dynamics is driven by a positive Levy process correlated with the 
jump process in the log-price of the asset, in the latter the volatility dynamics is 
governed by a time-changed Levy process. We shall concentrate in the present 
work on the second model we have just mentioned, the Bates model in which a 
Merton jump-diffusion model is combined with a stochastic volatility model of 
the Heston type. As R. Cont and P. Tankov have pointed out, the performance of 
the time changed Levy models in calibrating to market option prices are usually 
much better than those obtained in a BNS framework [I4j: in the last model, 
in fact, the possible implied volatility patterns are restricted by the requirement 
that the same parameter p characterize both the jumps in the returns and the 
volatility; on the other hand the calibration performances of the Bates model are 
comparable to those of the time changed Levy processes, "Thus the Bates model 
appears to be at the same time the simplest and the most flexible of the models" 
[17| . pag. 495. In the framework of option pricing via PDE (FIDE for models 
with jumps) several different approaches have been exploited both for stochastic 
volatility models and models with jumps. As far as finite element methods are 
concerned, we just quote the papers by Y. Achdou and N. Tchou [Ij and by N. 
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Hilber, A.-M. Matache and C. Schwab [14J for the first class of models and the 
papers by A.-M. Matache, T. von Petersdorff and C. Schwab [11] and by A.-M. 
Matache, P.-A. Nitsche and C. Schab [12] for the second. For models including 
both features the numerical methods proposed until now are much less. Some 
authors have considered finite-difference schemes for these models. D. Hilber, 
A.-M. Matache and Schwab have provided a finite-element approach to a large 
class of stochastic volatility models without jumps in [T3], including the Heston 
model. In the present paper we shall present a finite-element approach to option 
pricing in a Bates model framework. In the next section we'll recall the basic fea- 
tures of the Bates model, while in section 3 we'll provide the FIDE formulation 
of the option pricing in this model; in section 4 we shall present the variational 
formulation of the problem. In section 5 we'll describe the numerical method 
proposed and in section 6 we'll expose some comments on the results obtained. 
In section 7 we'll outline some conclusions and some futures perspectives of the 
present work. 



2 The Bates model 

In the Bates model the asset price evolution is given by: 

St = Soe''\ (1) 

where the log-returns X and its volatility Y satisfy the following stochastic 
differential equations: 

dXt = {a- Wt)dt + v^dW^/ + dZt, Xo = 0. (2) 

dYt = ^rj - Yt)dt + e^tdWl = y > 0, (3) 

where dW^ and dWf are two standard Wiener processes with correlation p 
and dZt is a Levy process of compound Foisson type . Let's assume for the 
parameters the following restrictions: 

OEM, -l<p<l, C>0, r/>0, 0>O (4) 

^2 < 2ir]. (5) 

The last requirement is in order to assure that the volatility process Y will 
never hit zero. We'll assume moreover E[Zl] < oo, this implying that the Levy- 
Khinchine representation formula for the process Z will be of the following type: 

k{z) =Cz + J (e^^^ - 1 - zx)U{dx) (6) 

where C = E[Zi] and U is the Levy measure of Z. We'll denote by fi{dx, dt) 
the jump measure of Z and by v{dt^ dx) its predictable compensator. We'll have 
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moreover iy{dx,dt) = U{dx)dt. We can write then: 

ft f + OO 



Ct + 



x{fi — u){dx, ds). 



10 J-oo 

The stochastic differential equation for the price S will be then: 

dSt = (a + nil))St^dt + St-VYtdWl+ 



St. 



— u){dx, dt). 



(7) 



(8) 



Remark 1 In the original model of 15], the process Z is a compound Poisson 
process, 

Nt 

Zt = ^ Ji, 



(9) 



i=l 



where N is a standard Poisson process with intensity A > and (Jj)j>i are iid 
N{'y, 6'^), with 7 = ln(l + k)— 6'^ /2. The corresponding cumulant function is in 
that case 

k{z) = A(e^^+'^'^'/2 _ ^10^ 

Remark 2 If Z = then we obtain Heston's stochastic volatility model from fBjj. 
If 9 = and y = r] we obtain Merton's jump diffusion model from 113^ with 
normal jumps. Consequently we might consider the Bates model as an extension 
of a Merton jump- diffusion model with stochastic volatility, or as an extension 
of a Heston volatility model with jumps in the returns. 



Lemma 1 The dynamics of the asset price process is given by 
dSt = (a + nil))St-dt + St^VYtdWl+ 

5t_(e^ - l){n-iy){dx,dt). 



In particular if 



a + k(1) 



(11) 



(12) 



the process S is a local martingale. 



Proof: This follows immediately from Ito's formula for general semimartin- 
gales. □ 

As in other affine stochastic volatility models with and without jumps, it 
is possible to obtain the characteristic function of the log-price in closed form. 
This characteristic function has been calculated by D. Bates in [5]; a detailed 
computation is provided also in [T7]; it is given by the following expression: 

"^^^g-(52M2/2+i(ln{l+m)-5V2)« _ 



exp 



cosh — + 



ipOu et 

smh — 

2 



exp 



(u^ + iu)vQ 



ecoth y + ^ — ipOu 



(13) 
(14) 
(15) 
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where: 

e = y^e^u^ + iu) + {^-ip9u)^ (16) 

Once the characteristic function of the log-price process is known in a closed 
form, the valuation problem for vanilla options can be easily solved by an FFT- 
related technique like that provided in the paper by P. Carr and D. Madan 

m- 

A quadratic approach to option hedhing in the Bates model has been sug- 
gested in [To]. 

3 Option pricing via PIDE approach 

By applying Ito lemma, and introducing the market price of risk tt, we obtain 
the following partial integro-differential equation for the price of a European call 
option C{S,y,t) in the framework of the Bates model : 

i^(v-y)-M^ + f'y^ + pOyS^+ (ir) 

/ + 00 
[C(5e^, y, t) - C{S, y, t)] W{dx) - rC = 
-oo 

with the following final condition at t = T: 

C{ST,yT,T) = max[ST-K,0] (18) 

and the following boundary conditions both in S and y: 

dC 

C(0,y,t) = 0,— (oo,y,t) = l (19) 
C7(5,oo,i) = 5, (20) 



oo 



C(5, 0, t) = V e-^*^CB5(S, t, K, an,fn), (21) 



n! 

n=0 



where Cbs{S, t, K, an,rn) ai'e the Black-Scholes values of the call options at time 
t for underlying price S and strike K with parameters 

an = n-i^jt, (22) 

f„ = r + A(l - e^+'^'/2^ + n(7 + 5V2)/i (23) 
By assuming lognormal jumps we have 

k(\) = A (e^-''/' - l) , (24) 
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where 7, (5^ are the expected value and variance respectively of the normal dis- 
tribution describing the jiuTips' size. 

The variables S, Y and t can assume values in the following domains: s G 
[0, +00), t G [0, +00) and y £ [0, +00) . 

Remark 3 The market price of risk tt can be obtained in different ways in 
the frame of general equilibrium models; consumption-based models give a risk 
premium proportional tu y. In the following, we'll assume without any lost of 
generality, that the market price of risk associated to the volatility is zero. The 
method can be generalized to a different choice in a straightforward way. 

Remark 4 By the commonly used substitution St = exp(Xt) , F{Xt,yt,t) = 
C{e^*,yt,t) the previous equation becomes: 

dF , y^dF ,dF 

-y^ + -9^y^ + pey-^-^+ (25) 
2 dx^ 2 dy'^ dydx 

/+00 
[F{X + u), y, t) - F{X, y, t)] W{du) -rF = 
-00 

W{dx) is the Levy density of the jumps. In the case of gaussian jumps for 
X (i.e. lognormal for S) it will be of the following form: 

W{du) = X-j=-^ eM^—^)du (26) 

It will be then a gaussian density with expected value 7 and variance 6^ multi- 
plied by the intensity A of the Poisson process. 

Remark 5 The boundary conditions for the new unknown F{x,y,t) are now 
the following. For X G : 

(27) 
(28) 



F{-oo,y,t) 


= 0, 


F(+oo,j/,t) 


= e", 


F{x,0,t) 


00 




n=0 


F{x,oo,t) 




The final condition is now: 





-At (At) 



CBs{e'',t,K,an,rn), (29) 

(30) 



F{Xt, yr, T) = max [e^^ - K, O] (31) 
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4 Variational Formulation 



The integro-differential equation given before can be written in the fohowing 
"divergence form": 



DF 

15t 

f +00 



+ V • [KVF] • 



/+00 
[F{X + u),y,t)-F{X,y, t)] W{du) -rF = 
-00 



(32) 



D d 

where the symbol -— = — — h (a • V) denotes the total derivative, the vector a 
^ Dt dt ^ ^ 

is given by: 

r-K{l) -y/2-p0/2 

i{ri-y)-e^/2 



a := 



(33) 



K 



(34) 



and the matrix K by: 

y/2 pey/2 
pey/2 e^y/2 

A variational formulation for the PDE arising in the Heston model has been given 
in pL5j, while the variational formulation for the PIDE in an Exponential Levy 
framework have been provided in [11], [12], [21], where existence and uniqueness 
of the solution of the variational problem associated with the differential and 
integro-differential equations respectively have been proved in suitable weighted 
Sobolev spaces and detailed analyses of both localization and discretization er- 
rors have been provided. Without performing the same analysis here we are 
going to proceed with the variational formulation for the present problem fol- 
lowing the same line. 

Introducing the following bilinear form: 



I 

Jn 



KVn • 'SJvdxdy 



bj{u,v) :-- 

u ( / [u{x + z) — u{x)] W{dz) \ dxdy 



(35) 



(36) 



it is possible to define a suitable discretization for our integro-differential 
problem [25l 



Remark 6 The equation 11 has degenerate coefficients both in the S and in 
the y variables; while the substitution x = log 5 removes the degeneracy in the 
S variable, this is still present for the y variable. In particular we want to 
discuss the boundary y = 0. If no restrictions on the model's parameters would 
be present, a condition on this boundary should be imposed in order to have a 
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well-posed problem and the correct condition is given hy\2^ . On the other hand, 
if the restriction on the parameters of the model given in section 2 holds, the 
variable y never hits that boundary. A closer inspection to the bilinear form 
associated to the problem suggests that when Green's lemma is applied to the 
l.h.s. of equation\2^ an "inflow" condition for our backward parabolic equation 
appears if the restrictions on the parameters hold and this in turn implies that 
the condition on y = need not to be imposed in order to have a well-posed 
problem. 



5 Numerical Approach 

In this section we will introduce a Finite Element Discretization of the above 
described PIDE. A similar approach for the Merton's and Kou's model has been 
introduced in [2]. 

5.1 Temporal Discretization 

The integro-differential equation can be written in the following form 
DF 

— +V-(KVF) + 

[F{x + u, y, t) - F{x, y, t)]W{du) - rF = 

-oo 

D d 

where the symbol -— = — + (a • V) denotes the total derivative. 

^ Dt dt ^ ' 

The time interval [0, T] will be discretized using a time step At such that 
t"+-^ = t^+At; moreover we will use the following notation = Fif^). Starting 
from this formulation we can obtain the following temporal discretization 

^ + V • (KVF"+^) - rF'^+V 

(38) 

(x + n, y, t) - F"+^(x, y, t)]W(du) = 

where F"(X) is the value of the price evaluated at the foot of the characteristic 
line at time and X is the solution of the following initial value problem 

^^^^^jAi^=a(r;X(r;t,x)) re(0,t) (39) 
X(t;t,x)=x. (40) 

This last ODE can be solved using either the implicit Euler method or a more 
accurate Runge-Kutta method. 

The characteristic Galerkin method, described above, is an Eulerian-Lagrangian 
approach and it is stable with a mild stability criterion allowing therefore the 
use of a large time step when appropriate. 
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5.2 Localization 

The infinite domain has to be reduced to a finite one: in particular we will 
consider a rectangular domain ft = [0,Xmax] x [0)2/moa;]- Using this reduced 
domain the boundary conditions have to be modifed accordingly in the following 
way: 

FiO,y,t)=0, (41) 
F{x, 0,t) = Y^ e-^'^FBsie"", t, K, an, a^), (42) 

n=0 

F{xmax,y,t)=e'^-\ (43) 
F{x,ymax,t)=e'. (44) 

Moreover the extrema of the integral term have to be reduced to finite values; 
to this aim we can use 



Ldown = - V \og{edV2^) - |7|, (45) 

and 

L„p = ^J-25■'\og{e5^/2^) + |7|. (46) 
5.3 Spatial Discretization 

The domain ft will be discretized using an unstructured triangular mesh. 
The discrete weak formulation reads as follows: find F^'^^ G such that 



/ F^"+Vdx- f F^"(X)V'dx + At6D(F^"+\V) 

Atbj{F;^+\ip) -Atf rF^+Vdx = VV' G 
Jn 



(47) 



where Vh is suitable functional space. In the present paper the unknown F will 
be approximated using Pi (linear) finite element i. e. 



NN 



FJi{x,y) = Y,FrUx,y), (48) 

i=l 

where A^A^ is the number of nodes of the triangulation and '^i{x,y) G Pi(r2). 
The problem has been solved using FreeFEM++. 

6 Numerical Results 

In this section we shall provide some comments on the numerical results ob- 
tained in order to assess the effectiveness of the proposed numerical method. In 
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Figure 1: Implicit Volatility Surface for parameter set SI. 

particular we'll present the implicit volatility surfaces for the following sets of 
parameters taken from [16j, where a suitable calibration methodology has been 
developed for a large class of stochastic volatility models with jumps. Moreover 
we'll provide a graphical comparison between the solution obtained with our 
finite element approach and that obtained by the usual method proposed by P. 
Carr and D. Madan based on the Fast Fourier Transform 

The range of the strike prices for the european call option considered is 
80 < K < 120 for Fig.l, 3, 4, while for Fig. 2 is 80 < K < 100. The maturities 
are between and 3 years. The initial value of the underlying S has been set 
S = 100. 



Set 




rj 


9 


P 


Si 


0.21568 


0.04937 


0.23828 


-0.44793 


S2 


0.33502 


0.033582 


0.26969 


-0.42404 


S3 


0.13279 


0.18193 


0.37518 


-0.59722 


5*4 


0.48443 


0.022097 


0.21903 


-0.40066 


Set 


k 


6 


A 




Si 


-0.11889 


0.17189 


0.13674 




S2 


-0.077973 


0.11048 


0.33785 




S3 


0.080396 


0.057373 


0.05218 




Si 


-0.12938 


0.16878 


0.15977 





All the volatility surfaces obtained exhibit both smiles and skews for short 
maturities. While the skew is quite pronounced for the set of parameters denoted 
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Figure 2: Implicit Volatility Surface for parameter set S2 




Figure 3: Implicit Volatility Surface for parameter set S3 
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0.6 



0.4- 




80 

Figure 4: Implicit Volatility Surface for parameter set S4. 

by 'S'2, 5*3 it seems less evident for the other sets Si, S4. The smiles appear less 
and less pronounced for longer maturities for all sets. The explanation of the 
different behavior exhibited by the volatility surfaces seems to be related more 
to the different values of the diffusion coefficient 9 of the volatility dynamics 
than to the other parameters, while the "leverage" coefficient p seems to be 
responsible of both the skews and the smiles appearing in the volatility surfaces. 
The behavior of the smiles in correspondence to the "at the money values" of 
the call option looks moreover quite realistic. 

The next two figures represent the solution obtained with the present fi- 
nite element method versus the solution obtained via a Fast Fourier transform 
method. While the latter is represented by the continuous line, the former is 
indicated by the dots. The set of parameters characterizing the model are those 
denoted before by 5i, ^4 respectively. The range of prices is 80 < 5 < 120. The 
calculation has been performed for T = 1, K = 100. 

A very good agreement between the solutions can be immediately recognized, 
although it looks quite evident that the solution obtained via the finite element 
method slightly overprices the call option for higher values of the underlying at 
time 0. 

In order to check the robustness of the present method with respect to pa- 
rameters changed, several trials have been performed corresponding to other 
sets of parameters belonging to an enlarged range of parameters and the results 
obtained look very close to those just presented. 

The CPU time required for a calculation of the call option price for a single 
value of the underlying is about 3 minutes, while that required for the volatility 
surfaces shown here is about 3 hours on a 2Ghz Centrino Duo with 1MB RAM. 
The 2D mesh used in the computation is composed by 2634 triangles and 1398 
nodes. 
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7 Concluding Remarks 



We have presented a finite element approach to the european option valuation 
problem formulated via a partial integro-differential equation. Several choices 
of the functional space suitable for the spatial discretization are possible and 
we have made the most simple choice in order to obtain a fast and efficient 
algorithm. Other choices have been made by some authors in similar contexts, 
like in [S], J2]) where Wavelets have been used in an extensive way to produce 
an accurate algorithm, but the convergence of the method with this choice seems 
to be not very fast. In the present context, with a single underlying, the choice 
of piece-wise polynomial function spaces seems to perform slightly better. 

A natural development of the present analysis will be the variational formu- 
lation and the implementation of a finite element method for American option 
pricing with the Bates model. This problem can be formulated as a free bound- 
ary problem for the same Partial Integro-Differential equation given before. The 
finite element approach seems to be the most promising way to obtain fast and 
accurate algorithms for this problem. This will be the subject of future investi- 
gations. 
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